function res = adp(f,a,b,e)
% ADP - 基于Simpson公式的自适应积分法求函数f从a到b的积分。
%       e为最大容许误差
    if nargin < 4
        e = 3*eps;
    end

    if e<3*eps
        error('Epsilon too small')
    end
    res = adpt(f,a,b,e);
end

function res = adpt(f,a,b,e)
% ADPT - 基于Simpson公式的自适应积分法求函数f从a到b的积分。
%       e为最大容许误差
    
    H = b-a;
    h = H/2;
    m = a+h;
    l = (a+m)/2;
    r = (m+b)/2;

    inte1 = H*(f(a)+f(b)+4*f(m))/6;
    
    inte2 = h*(f(a)+f(m)+4*f(l))/6+h*(f(m)+f(b)+4*f(r))/6;
    
    E = abs(inte1-inte2);
    
    if E < e
        res = inte2;
    else
        res = adpt(f,a,m,e/2)+adpt(f,m,b,e/2);
    end    
end

